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Abstract 

We study front speeds of curvature and strain G-equations aris- 
ing in turbulent combustion. These G-equations are Hamilton- Jacobi 
type level set partial differential equations (PDEs) with non-coercive 
Hamiltonians and degenerate nonlinear second order diffusion. The 
Hamiltonian of strain G-equation is also non-convex. Numerical com- 
putation is performed based on monotone discretization and weighted 
essentially nonoscillatory (WENO) approximation of transformed G- 
equations on a fixed periodic domain. The advection field in the com- 
putation is a two dimensional Hamiltonian flow consisting of a peri- 
odic array of counter-rotating vortices, or cellular flows. Depending on 
whether the evolution is predominantly in the hyperbolic or parabolic 
regimes, suitable explicit and semi-implicit time stepping methods are 
chosen. The turbulent flame speeds are computed as the linear growth 
rates of large time solutions. A new nonlinear parabolic PDE is pro- 
posed for the reinitialization of level set functions to prevent piling up 
of multiple bundles of level sets on the periodic domain. We found that 
the turbulent flame speed st of the curvature G-equation is enhanced 
as the intensity A of cellular flows increases, at a rate between those of 
the inviscid and viscous G-equations. The st of the strain G-equation 
increases in small A, decreases in larger A, then drops down to zero at 
a large enough but finite value A* . The flame front ceases to propagate 
at this critical intensity A*, and is quenched by the cellular flow. 
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1 Introduction 



Front propagation in turbulent combustion is a nonlinear and multiscale dy- 
namical process [351 113 EH E31 E3 El EE ES] • The first principle based ap- 
proach requires a system of reaction-diffusion-advection equations coupled 
with the Navier-Stokes equations. Simplified models, such as the advec- 
tive Hamilton-Jacobi equations (HJ) and passive scalar reaction-diffusion- 
advection equations (RDA), are often more efficient in improving our under- 
standing of such complex phenomena. Progress is well documented in books 
[351 EHl [37] and research papers PSITJIHlIIOlIIlISlEZlEIlESlIMlEail^ 
among others. 

A sound phenomenological approach in turbulent combustion is the level 
set formulation [27] of flame front motion laws with the front width ignored 
|29j . The simplest motion law is that the normal velocity of the front (V n ) 
is equal to a constant sl (the laminar speed) plus the projection of fluid 
velocity V(x, t) along the normal n. The laminar speed is the flame speed 
due to chemistry (reaction-diffusion) when fluid is at rest. As the fluid is in 
motion, the flame front will be wrinkled by the fluid velocity. Under suitable 
conditions, the front location eventually moves to leading order at a well- 
defined steady speed st in each specified direction, which is the so-called 
"turbulent burning velocity" [29J. The study of existence and properties of 
turbulent flame speed st is a fundamental problem in turbulent combustion 
theory and experiments [65\ [3T| [2"§] . Let the flame front be the zero level set 
of a function G(x,t), then the normal direction is DG/\DG\ and the normal 
velocity is —Gt/\DG\. (D: spatial gradient.) The motion law becomes the 
so-called G-equation in turbulent combustion [351 ESj : 

G t + V(x,t)-DG + s L \DG\ =0. (1.1) 

Chemical kinetics and diffusion rates are all included in the laminar speed 
sl which is provided by a modeler. Formally under the G-equation model, 
for a specified unit direction P, 

s T (P) = - lim (1.2) 



Here G(x,t) is the solution of equation (1.1) with initial data G(x, 0) = 
P ■ x. The existence of st has been rigorously established in |38j and [5] 
independently for incompressible periodic flows, and [22j for two dimensional 
incompressible random flows. 
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As fluid turbulence is known to cause stretching and corrugation of 
flames, additional modeling terms may be incorporated into the basic G- 



equation (1.1). In this paper, we shall study turbulent burning velocity st 
of such extended G-equation models involving strain and curvature effects. 
The curvature G-equation is: 

G t + V(x, t)-DG + s L \DG\ = ds L \DG\div . (Gc) 

which comes from adding mean curvature term to the basic motion law. 
The curvature dependent motion is well-known, see \28\ [27] and references 
therein. If the curvature term is further linearized [IT], we arrive at the 
viscous G-equation: 

G t + V{x,t)-DG + s L \DG\=ds L AG, (Gv) 

which is also a model for understanding numerical diffusion [27] . The strain 
G-equation is: 

Gt + V(x, t)-DG+ (a L + d DG '%Q 2 DG > ) \DG\ = ds L \DG\div . 

(Gs) 

The strain term n ■ DV ■ n will be derived and analyzed later. The formula 



(1.2) formally extends to (Gc) and (Gs). A complete mathematical theory 



of their existence is lacking at the moment. Helpful empirical observations 
from experiments [H [31] are: (i) When the flame front is wrinkled by the 
advection, the interface area increases and st increases (called "enhance- 
ment"), (ii) However, turbulent flame speed cannot increase without limit, 
and the growth rate may be sublinear in the large intensity limit of the ad- 
vection (called "bending"), (iii) When the advection is strong up to certain 
level, the reactant totally scatters. The reaction then fails and the flame 
front extinguishes (called "quenching"). 

We aim to understand and quantify these nonlinear phenomena in the 
context of curvature and strain G-equations and cellular flows where st is 
related to the corrector (cell) problem of homogenization theory for which 
several mathematical results are available. The cellular flow is a two dimen- 
sional incompressible flow: 

V = V^U = (-H y ,n x ) , U = ^- sin(27rx)sin(27ry), (1.3) 

where A is the amplitude of the flow. By parameterizing st as a function 
of A, we are interested in the behavior of st as A increases in G-equations 
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( 1 . 1 ) , ( Gc ) , ( Gv ) , ( Gs ) . The streamlines of the cellular flow consist of a pe- 



riodic array of hyperbolic (saddle points, separatrices) and elliptic (vorti- 



cal) regions. For inviscid G-equation (1.1), it is known |26| [TJ [23] that 
st = 0(A/ log(A)), where the logarithmic factor is due to slow-down of 
transport near saddle points. For viscous G-equation, we recently proved 
[16] that st = 0(1) as A S> 1 at any fixed positive viscosity (d > 0). The 
dramatic slowdown (strong bending) is due to the smoothing of the level 
set function G by viscosity, and the uniform bound of ||-DG|| L i . Less is 

loc 

known about the growth rate of st for curvature and strain G-equations. 
The curvature term only provides partial smoothing, hence the slowdown 
(bending) is weaker in general than the regular smoothing by viscosity. For 
shear flows, we showed [16] that the linear growth rate lim J 4_ >00 st/A is same 
as that of the inviscid G-equation. The effect of strain term is more difficult 
to analyze, as it is highly nonlinear in G and can take both signs. It also 



changes the type of the Hamiltonian of G-equation from convex in (1.1) to 



non-convex in (Gs). For shear flows, the strain term always slows down st 
[39]. 

We shall first approximate the G-equations by a monotone discrete sys- 
tem, then apply high resolution numerical methods such as WENO (weighted 
essentially non-oscillatory finite difference methods |131 [2"7] ) with a combi- 
nation of explicit and semi-implicit time stepping strategies, depending on 
the size and property of dissipation in the equations. The computation is 
done on transformed G-equations over a periodic domain to avoid the need 
of excessively large computational domains to contain potentially fast mov- 
ing fronts. We also devise a new reinitialization equation on the periodic 
domain to prevent the level sets from piling up during time evolution. A 
nonlinear diffusion term is added to the standard reinitialization equation 
(Chapter 7, [27]) to perform reinitialization on multiple bundles of level sets 
often encountered during long time computation. An iterative method of 



computing st of the viscous G-equation (Gv) works well based on the cor- 



rector equation of homogenization, if the viscosity d is above a certain level. 



Our main findings are: (1) The curvature G-equation (Gc) always en- 
hances st as A increases; the amount of enhancement is smaller than that 
of the inviscid G-equation (1.1), larger than that of the viscous G-equation 
(Gv). For small enough d, the st of (Gc) behaves similarly to that of the 
inviscid G-equation (1.1), or weak speed bending. For large enough d, the 
st of (Gc) behaves similarly to that of the viscous G-equation (Gv), or 
strong speed bending. (2) The st is a monotone decreasing function of d for 
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both curvature and strain G-equations, (Gc) and (Gs). (3) For the strain 
G-equation (Gs) with fixed d > 0, st first increases with A, then decreases 
in A, and drops down to zero at finite A (front quenching occurs). 



The paper is organized as follows. In section 2, we give a brief deriva- 
tion of G-equation models and an an overview of analytical results of the 
turbulent flame speeds. In section 3, we introduce numerical scheme for 
each G-equation. We also discuss how to perform reinitialization in periodic 
domain. In section 4, we present and interpret the numerical results. Con- 
cluding remarks are in section 5. In the two appendices, we show a formula 
of surface stretch rate in advection and a convergent iteration scheme of st 
based on the corrector problem of homogenization. 

The work was partially supported by NSF grants DMS-0911277 (JX) 
and DMS-0901460 (YY). YL thanks the Department of Mathematics of UC 
Irvine for a graduate fellowship. The authors would like to thank Professor 
John Lowengrub for helpful conversations on interface computations. 



2 Derivation and Analysis of G-equations 
2.1 G-equations 

In the thin reaction zone regime and the corrugated flamelet regime of pre- 
mixed turbulent combustion (pp. 91-107, Chapter 2, [29]), the flame front is 
modeled by a level set function: {(x, t) : G(x, t) = 0}, which is the interface 
between the burned area {G < 0} and the unburned area {G > 0}. See 
[27j for an introduction on level set methods in a broad context. The unit 
normal direction is n = DG/\DG\ and the normal velocity is — Gt/\DG\. 
(D: spatial gradient.) The simplest motion law is that the normal velocity 
of the interface is the sum of a constant sl (called laminar flame speed) and 
the projection of fluid velocity V(x, t) along the normal direction. The sl is 
well-defined if the reaction zone is much larger than the smallest turbulent 
length scale (the Kolmogorov scale), as in the corrugated flamelet regime 



[29] . In terms of G, the law is the so-called G-equation (1.1). A linear ver- 
sion dated back to [19]. The trajectory of a particle x(t) on the interface 
satisfies: 

= V(x,t) + s L n. (2.1) 

The G-equation or level set framework is a popular and robust phe- 
nomenological approach. The motion law is in the hands of a modeler based 
on theory and experiments. Various nonlinear effects may be built into the 



5 



basic model (1.1). For example, turbulence is known to cause stretching of 
flame fronts. It was shown in [21\ [2U] that the flame stretch rate may be 
added as a first order correction term on the laminar flame speed: 

, 1 da , . 

s L = s L -d- — , (2.2) 
a at 

where a is the surface element area of the level set and d is called the 
Markstein diffusive number. If the flame stretch rate is positive, the reactant 
on the flame front scatters and the burning reaction slows down. By a 
kinematic calculation (see appendix A for details), the flame stretch rate is: 

1 da , 
-—=S + s l k, (2.3) 
a at 

S = —n ■ DV ■ n , k = div(n), 
where S is called the strain rate and n is the mean curvature of level set. 



Replacing sl by §l in ( |2.1[ ), we have the strain G-equation (Gs). In the 
thin reaction zone regime (section 2.6, pp 104-107 [29]), Kolmogorov scale 
eddies enter the reaction zone, and cause unsteady perturbations of laminar 
speed sl. The (sl — dS) term and the eddy effects are lumped together as a 
fluctuating quantity (denoted by sl, s m |29j ) which however is on the order 
of Si based on direct numerical simulation data. If we approximate sl s by 



sl and keep the curvature term, the curvature G-equation (Gc) follows. 



Remark 2.1 In previous works \41\ Wj, sl is modified to remain positive: 

f Ida \ ( d lda\ 

s L = raax< s L - d- — ,0 > , s L exp — . 



a dt J \ sl a dt J 

However, these modifications restrict the curvature or strain effect in the 
strong advection scheme. The bending or quenching effect may either weaken 
or disappear. 

2.2 Turbulent Burning Velocity 

We discuss how to evaluate turbulent flame speeds in G-equation models. 



For simplicity we consider the inviscid G-equation (1.1) only, and the for- 
mulation extends to other G-equations. 

Given a unit vector P G W 1 and suppose the flame front propagates 
in direction P. Let the initial flame front be {P ■ x = 0} and consider 
G-equation with planar initial condition: 

j G t + V(x, t)-DG + s L \DG\ =0 in M n x (0, oo) 

\ G(x,0)=P x onrx{t = 0} ' (2 '" 
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Assume V(x,t) is spatially periodic. If we write G(x,t) = P ■ x + u(x,t), 
then u(x, t) is also spatially periodic and solves the following periodic initial 
value problem: 

j u t + V(x, t)-(P + Du) + s L \P + Du\ =0 in T n x (0, oo) 

\u(x,0) = on T" x {t = 0} ' ^'^ 



Hence in numerical computation of (2.4) we can reduce the spatial domain 
from M. n to [0, l] n by imposing the affine periodic condition: 

( G t + V(x,t)-DG + s L \DG\ =0 in [0, If x (0, oo) 

\ G(x,0) = P-x on [0, 1]" x {t = 0} ' U '°' 

G{x + z,t) = G(x, t)+P-z, x G [0, l] n , 2 G Z n . (2.7) 

Now we focus on P = ei = (1, 0, • • • ,0), then G(x, t) = x\ + i) is 
periodic in X2, ■ ■ ■ ,x n . Consider the stripe domain M x [Ojl]™" 1 , and the 
burned area at time t is {x G R x [0, l] n_1 : G(x,t) < 0}. Denote A(t) 
the volume that burned area has invaded during time interval (0, t), then 
turbulent flame speed is the linear growth rate of A(t): 

Ait) 1 f 

ST = AToo ~T = tiToo t / r , 1 (X{GGM)«» - X{G(x,o)<o}) da;- (2-8) 

(x- indicator function.) Note that G(x, 0) = x\ and G(x+ei, t) = G(x, t)+l, 
then A(t) and hence st can be evaluated by G or u in [0, l] n : 

s T = lim — / [G(x,t)]dx= lim — / [xi + u(x, t)] dx. (2.9) 

t^+oo t 7[o,l]» t^+oo t J[0,1]" 

([•]: floor function.) In |41j the initial condition is chosen as G(x, 0) = <f>(x\) 
with 4> '■ R — ^ R a smeared-out signed function, and the computational 
domain is [a, b] x [0,1]. If the zero level set travels a long distance, the 
length of the domain (6 — a) needs to be large enough to contain the level 
set. To study a fast moving flame front and its long time behavior, the 
computational domain will be very large. Instead we choose G(x, 0) = x\ 
and reduce the computational domain to [0, 1] x [0, 1]. The st is the same 
from either initial data. 

Another way to find turbulent flame speed is via the framework of pe- 
riodic homogenization |14^ I12j . Assume V = V(x) be time-independent 
periodic flow and consider the so-called corrector problem: given any vector 
P G R n , find a number H (the effective Hamiltonian) such that the equation 

V{x) ■ (P + Du) + s L \P + Du\=H, x G T n (2.10) 
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has a periodic solution u(x). If (2.10) is solvable, then G-equation has the 
following stationary solution: 



G(x, t) = -Ht + P-x + u(x), 



(2.11) 



and H is exactly the turbulent flame speed. The corrector problem is well- 
posed for viscous G-equation [15] , and can be used to compute st iteratively 
when viscosity is not too small (see section 3.4 and Appendix B). However, 



(2.10) for inviscid G-equation may not have exact solutions due to lack of 



coercivity of G-equations, only approximate solutions exist (38]. It is also 
an open question in general whether it has solutions if the curvature or 
strain term is present. The more general and robust characterization of st 
is simply the linear growth rate of G or u at fixed x: 



s T 



lim 

t— >+oo 



-G(x,t) —u(x,t) 
= hm . 

t t->+oo t 



(2.12) 



which we shall adopt for curvature and strain G-equations in this paper. 
Indeed (2.9) and (2.12) are consistent when P = e±, but (2.12) can be 

See 



used for any direction P. See [30] for earlier work on computing effective 
Hamiltonian of coercive Hamilton-Jacobi equations along this line. 



3 Numerical Methods 

We discuss the numerical schemes for G-equations. We employ the Hamilton- 
Jacobi weighted essentially nonoscillatory (HJ WENO) scheme and the total 
variation diminishing Runge-Kutta (TVD RK) scheme in higher order spa- 
tial and time discretization respectively. See [131 [32j [27] for details of the 
schemes. 



3.1 Inviscid G-equation 



Inviscid G-equation (1.1) is a Hamilton-Jacobi equation with Hamiltonian 

H{p) = V(x,t)-p + s L \p\. (3.1) 



The forward Euler time discretization of (1.1) is 

Qn+1 



G n 



At 



+ H n (G x , G x , G v , G v ) — 0, 



(3.2) 



where H is the numerical Hamiltonian of (3.1 ) and G x (G x ) denotes the left 
(right) discretization of G x . For the stability of the numerical scheme, H = 
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H(p x ,p£ ,p y ,Py) is chosen to be consistent and monotone [9]. Here con- 
sistency means that H(p x ,p x ,p y ,p y ) = H(p x ,p y ), and monotonicity means 
that H is nondecreasing in p~ , p~ and nonincreasing in p x ,p y ■ 



Write V = (Vi,V 2 ) in (3.1): 



H(p x ,Py) = ( Vx + ) p x + I V 2 + S L 



Py 



Py 



When the velocity field dominates the normal velocity, upwinding direction is 
determined by the velocity field. For example, if V\ > sl, then V\ + slPx/\p\ 
is always positive and p x is approximated by p~ . However, if the velocity 
field and the normal velocity are comparable, it is hard to determine the 
upwinding direction. In this case we treat both terms separately: for the 
velocity field term, we apply upwinding scheme; for the normal velocity term, 
we apply Godunov scheme. Since both schemes are monotone, their sum is 
again monotone. In summary, we have the following monotone numerical 



Hamiltonian of (3.1): 



H{p-, P +, P -,P+) = V lP V x d + V 2P V / 1 + S Lx /(p^)2 + (pNorV 



where 



and 



p x , if Vi > 
Pt ,ifFi<0 



Nor\2 



Nor\2 



, if V 2 > 
,ifF 2 <0 

f ( Px ) 2 ,iSV 1 >a L 

max (max(p~, 0) 2 , min(p+, 0) 2 ) ,if | Vx] < s_l 
I (Pt? MVk-sl 

{ (Py? AtV 2 >S L 

max (max(p~, 0) 2 , min(p+, 0) 2 ) , if \V 2 \ < sl 

{ (pi) 2 mv 2 <s l 



(3.3) 



For the accuracy of the numerical scheme, we apply WEN05 scheme to 
approximate the spatial derivatives and RK3 scheme in forward Euler time 
discretization. The time step restriction (CFL condition) is 



At 



IK II +s L \\V 2 \\ +s L | 



Ax 



Ay 



(3.4) 



(|| • || : maximum norm.) Overall the scheme gives nearly fifth order spacial 
accuracy in smooth regions of solutions, and third order accuracy in time. 
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Remark 3.1 Compared with the standard schemes (LF, LLF, RF, etc. See 
chapter 5 of \27j ), our choice of numerical Hamilton is easy to implement, 
and no extra artificial diffusion is added to satisfy the monotonicity. 



3.2 Curvature G-equation 



In forward Euler scheme of curvature G-equation ( Gc ) , the curvature term 
in two dimensional space is 



|£>G|div 



DG \ GyG xx — 2G x GyG X y + G x G yy 



\DG\ J G 2 x + G 2 y 



and is discretized by central differencing |28j . Since central differencing 
gives only second order accuracy, we apply WEN03 scheme to evaluate the 
numerical Hamiltonian (3.3) and RK2 scheme in time step discretization. 
The time step restriction is 

A , (\m±^L , ra+sL M 2 SL d 2 SL d \ 

At {-^ x — + ^ y — + (A^ + (A y j) <1 - (3 - 5) 

When d is large (3> Ax), the time step size for forward Euler scheme is 
very small At = 0((Ax) 2 ). To alleviate the stringent time step restriction, 
we decompose the curvature term as follows: 

\DG\div f p§|) =AG- AooG (3.6) 

- ir j- r \ - ^^ xx + 2G x G y G xy + G 2 y G yy 

where Aqo is the infinity Laplacian operator. If we apply backward Eu- 
ler scheme on AG and forward Euler scheme on AooG, then we have the 
following semi- implicit time discretization scheme for (Gc): 



At 



+ V(x, t n ) ■ DG n + s L \DG n \ = ds L (AG n+L - A^G 71 ), (3.7) 



whose time step restriction is same as inviscid G-equation (3.4). Note that 
for implicit scheme each time step is more expensive. Hence if d is small 
(~ Ax), the forward Euler scheme is still the better choice. 

Another cause of small time step is when || V|| is large. However we can- 
not move the advection term into implicit scheme as in standard advection- 
diffusion equations. The curvature G-equation is essentially of hyperbolic 
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type rather than of parabolic type. Even involving second order derivatives, 
the curvature term is dissipative only along the tangential plane of the level 
set and so cannot stabilize the advection term. 



Remark 3.2 The curvature term and the infinity Laplacian operator in 
higher dimensional space are 



\DG\div 



\\DG\J 



''.I 



G Xi G x 

~\dgJ< 



Gt - Gt ■ 



3.3 Strain G-equation 



For strain G-equation (Gs), the Hamiltonian becomes 



H(p) = V(x,t)-p + (s L -dS)\p\ , S 



p • DV ■ p 



(3.8) 



If we apply upwinding scheme on V ■ p, then it suffices to find a monotone 
scheme for (sl — dS)\p\. First we approximate p to obtain S, and next we 
evaluate \p\ by Godunov scheme according to the sign of (s_l — dS). Then 
we obtain the following monotone numerical Hamiltonian of (3.8): 



H(p-, P i,Py,pt) = V lP l el + V2P™ + (S L - dS)J{p^r)2 + ( p iV r )2) (3.9) 



where p^ el , Py el are same as in ( 
5 with p evaluated by central differencing, and 

(PS 



3.3), S is the numerical approximation of 



max(max(p x , 0) 2 , min(p+, 0) 2 ), if (sl — dS) > 
max(min(p~, 0) 2 , max(p+, 0) 2 ), if (sl — dS) < 



max(max(p J/ , 0) 2 , min(p,j~, 0) 2 ), if (sl — dS) > 
max(min(p~, 0) 2 , max(p+, 0) 2 ), if (sl — dS) < 



Remark 3.3 For cellular flow (IS), the strain rate can be simplified as 

(Gy — G 2 ) 



-2irA cos(27ra;) cos(27ry) - 



\DG\< 



Then (sl — dS) is always positive if 2irAd < sl- 
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3.4 Viscous G-equation 



When d is small, viscous G-equation ( Gv ) is advection dominated and should 
be treated like a hyperbolic equation. Similar to curvature G-equation, for 
spatial discretization, we apply WEN03 scheme on numerical Hamiltonian 
(3.3) and central differencing on the diffusion term. For time step discretiza- 
tion, we apply RK2 forward Euler scheme. 

When d is large enough, we consider the following semi-implicit scheme: 

rin+l rin 

At + V(x, t n+1 ) ■ DG n+l + s L \DG n \ = ds L AG n+1 . (3.10) 

Here the advection and diffusion terms are discretized by central differenc- 
ing, and the normal direction term is discretized by Godunov and WEN03 
scheme. Since there is no time step restriction from both advection and 
diffusion terms, the time step constraint for (3.10) is 

\Ax AyJ 

When V = V{x) is periodic, mean zero and incompressible, turbulent 
flame speed may also be obtained from the corrector problem: 

- ds L Au + V(x) ■ (P + Du) + s L \P + Du\ = H, x £ T n , (3.11) 

which has a unique (up to a constant) classical solution and H = sl J t « + 
Du\dx. When d is large enough, the following iteration scheme converges: 

-ds L Au (k+1) + V(x)-Du (k+1) = - s L \P + Du {k) \ - V(x)-P, x G T n , 

H^ = s L [ \P + Du {k \dx. (3.12) 

A convergence proof is in Appendix B. To solve (3.12) numerically as an 
elliptic equation, all operators are discretized by central differencing. 



3.5 Reinitialization 

When the flame front travels very fast, the level set function becomes very 
flat. When the motion of the flame front nearly stops, the level set function 
becomes very sharp. In either case the computational error will increase, 
and the level set may not be well captured. Hence reinitialization needs to 
be applied regularly to keep the level set function approximately equal to 
the signed distance function near the level set. 
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(e) 
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1 1 





Figure 1: (a) Graphs of S(4>) and D((f>). (b) Contour plot of the testing 
4>(x,t). (c) Reinitialized 4>(x,t) without smoothing at t = 0.2. (d) Reini- 
tialized (f)(x, t) at t = 0.2 with one smoothing iteration every time step, (c) 
Reinitialized <j>{x, t) at t = 0.2 with 10 smoothing iterations every time step. 
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The standard reinitialization equation is 



4> t + Sm\D^\ -1) = , S(0) = sgn( 



(3.13) 



which spreads out the signed distance from the level set {cj>(x,t) = 0}. The 
function S : R — > R can be mollified to improve the numerical accuracy, see 
chapter 7 of [27] for details. 

ei, 4>(x,t) must satisfy 
Figure [I] (b) for 



To perform reinitialization on (2.6) with P 

(x + ei,t) : 



(x, i) + 1 and be periodic in X2, • 
an example of 4>(x, t) using contour plot. To maintain the spatial periodicity, 



we modify (3.13) as follows: 



<f> t + S(<t>)(\D<j>\-l) = 0, 



(3.14) 



where S is a 1-periodic function and S((f>) = sgn(<^>) for (f> G [—1/2,1/2]. 
See Figure [lj (a) for the graph of the mollified version of S((f>). In numerical 
computation, ( |3.14[ ) is discretized by WEN05 scheme and RK3 schemes 
with time step At = Ax. 

However, Figure [TJc) shows that (3.14) spreads out distances from both 
{4> = 0} and {eft = 1}. As a result, <fi(x,t) is squeezed near {<f> = 1/2}. The 
computation grinds to a halt when <f)(x, t) becomes too sharp. To avoid this 
problem, we consider the following nonlinear diffusion equation: 



4> t = cD(<f>)A<f>, 



where c is some positive constant and D : R 
satisfying D{4>) = for € [— e, e] and D{4>) 
Figure [l](a) for the graph of D{4>). Equation (3.15) smooths 



(3.15) 

1 is a 1-periodic function 
l_for 4> G [2c, 1 - 2e]. See 
x, t) in the 



region away from the level set. In summary, we combine (3.14) and (3.15) to 



obtain the following reinitialization equation for the transformed G-equation 
d2i6l) with P = ei: 



<j> t + S{(j>){\D ( j>\-l) = cD{(j>)A4>. 



(3.16) 



In actual computation, we do not solve (3.16) accurately because the dif- 
fusion term reduces the time step to At = 0((Ax) 2 ). Instead, we alternate 



between (3.15) and (3.14). Approximate (3.15) by the simple iteration: 



H+l,j + 4>i-l,j + 4>i,j+l + 4>i,j-l) 



. (3.17) 



The iteration (3.17) is repeated a few times in each time step of the numerical 



scheme of (3.14). This way, the time step remains At = O(Ax). See Figure 
[T|d) and (e) for an illustration of the smoothing effect. 
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4 Numerical Results 



We consider all G-equations (l.l),(Gc),(Gv),(Gs) in two spatial dimensions 



with P = e\ and sl = 1- The velocity field V(x,t) is chosen to be cellular 



flow (1.3) with various values of the intensity A to study the growth rate of 
turbulent flame speed. Also the Markstein number d is varied to study the 
curvature and strain effect. 



First we solve the periodic initial value problem (2.6) for G(x,t) on 
[0, l] 2 . Then by (2.7) we construct the solution G(x, t) in some stripe domain 



[a, b] x [0,1] and obtain the level set {G(x,t) = 0}. The computational 
domain is [0, 1] x [0, 1] with grid points up to 400 x 400. 

Figure [2] shows the graphs of G(x,t) for inviscid, curvature, and viscous 
G-equations at t = 1 with A = 4,8,16 and d = 0.1. When A is large, 
the graph of G(x, t) has cone shape in each cell. Due to the curvature 
effect, G(x, t) is less irregular and the cone formation is slower. The regular 
viscosity makes G(x,t) even smoother. 

Figure [3] shows the contour plots of G(x,t) for inviscid and curvature G- 
equations. When the level set merges, shock waves occur and the derivative 
of G(x, t) is discontinuous across the shock wave. We observe that the shock 
wave is of spiral shape in each cell, especially at d = 0.1, A = 16. 

Figure [4] shows the propagation of the flame front for inviscid and cur- 
vature G-equations at A = 32 and d = 0.1. When A is large, the flame front 
of the inviscid G-equation travels faster along the boundaries of the cells 
with bubbles formed behind. The flame front spirals inside the cells, and 
the bubbles shrink in the wake. If the curvature effect is added, the flame 
front is concave when traveling along the boundaries. The curvature term 
slows down front propagation yet the wake bubbles shrink faster. 

Figure [5] shows the time derivative function of A(t) and G(x = 0,t) for 
inviscid, curvature and viscous G-equation with ^4 = 8 and d = 0.1. After 
a short time interval, A'(t) behaves like a periodic function. Hence we can 
approximate st by taking the average of A!{t) over a periodic time interval: 



T2 — T\ J Tl T2 — T\ 

See Figure [5] for examples of selections of T\ , T2 . So we don't need to use 



( 2.9 ) and perform large time simulation in order to approximate st correctly. 

Next we consider the behavior of G'(x,t) in time for fixed x. (': d/dt.) 
For inviscid G-equation, G'(x, t) behaves like a periodic function after a short 
time, hence we can evaluate st by the same method as above rather than 



using (2.12). For viscous G-equation, the dissipation term causes damping 
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(a) Inviscid G-equation 




o 



(b) Curvature G-equation 




p 



(c) Viscous G-equation 




DO 



Figure 2: Graphs of G(x,t) at t = 1 for inviscid, curvature, and viscous 
G-equations in cellular flow with A = 4,8, 16 (left to right) and d = 0.1. 
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(a) Inviscid G-equation 




(b) Curvature G-equation 





Figure 3: Contour plots of G(x, t) at t = 1 for inviscid and curvature G- 
equations in cellular flow with A = 4,8, 16 (from left to right) and d = 0.1. 

(b) Curvature G-equation 



Inviscid G-equation 


I : 







m 



-•mm 


• 
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• • 





Figure 4: Propagation of flame front in time for inviscid and curvature G- 
equations with A = 32, d = 0.1 and t = 0.1, 0.2, 0.3, 0.4, 0.5. 
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(a) Inviscid G-equation (b) Curvature G-equation (c) Viscous G-equation 
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Figure 5: Plots of A'(t) and G'(0,t) for inviscid, curvature and viscous G- 
equation with A = 8 and d = 0.1. 



in G'(x,t). Hence G'(x,t) converges to — st in time, and G(x,t) converges 



to the stationary solution (2.11 ). For curvature G-equation, however, we see 



only slight damping in G'(x,t). 

Figure [6] shows function G'(0, t) with different grid sizes. For inviscid 
G-equation, the numerical scheme is higher order accurate, and the artificial 
dissipation is well minimized. Hence damping is hardly observed even on 
coarse grid. For curvature G-equation, the numerical scheme is second order 
accurate, and the curvature term may be incorrectly evaluated at shock 
wave. Hence damping effect is very strong on coarse grid, and we must use 
fine grid to reduce the artificial diffusion. 

We denote s l ™ , s™ r , s^ r the turbulent flame speeds for inviscid, 
curvature, viscous, strain G-equations respectively. We also denote them as 
functions of either the flow intensity (^4) or the Markstein number (d). Note 

s™ r (d) = s^ s {d) = s^ r (d) = sl, and when 



we have s l ™ 



that when A 

d = 0we have sp v (A) = s^ r (A) = s^ s (A) = s% r (A). 

Figure Q a) shows the graphs of s™ v {A), s™ r {A) and s1f s (A) with d = 
0.1. The numerical results indicate that they all increase as A increases and 

sT{A) < sf r {A) < sip* (A). 

Figure Qb) shows the graphs of sf v (A) and sf T {A) with d = 0.1,0.2,1. 
We used the forward Euler scheme for d = 0.1 and semi-implicit scheme for 



18 



(a) Inviscid G-equation (b) Curvature G-equation 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.3 1 

* i 

Figure 6: Plots of G"(0, t) for inviscid and curvature G-equation with A = 8, 
d = 0.1 and grid sizes 100, 200, 400. 




Figure 7: (a) Plots of st = st(A) for inviscid, curvature and viscous G- 
equations with d = 0.1. (b) Plots of st = st(A) for curvature G-equation 
with d = 0.1, 0.2 and 1. 
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Strain G-equation 


(d = 0.01) 


(b) Strain G-equation (d 
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Figure 8: Propagation of flame front in time for strain G-equation with 
A = 32, d = 0.01, 0.02 and t = 0.3, 0.6, 0.9, 1.2, 1.5. 



d = 0.2 and 1. It is known that sf v (A) = 0(A/ log A) and s^ is (A) = 0(1). 
However, the precise asymptotic behavior of Sj? r (A) as A — > oo remains 
open. The growth scaling of s™ r (A) is not conclusive from the range of A 
we simulated. 

Figure [8] shows the propagation of the flame front for strain G-equation 
with A = 32 and d = 0.01, 0.02. Near the corner of the cell, the velocity 
field is weak (|V(a;)| ~ 0) yet the strain rate is strong (|«S| ~ 2irA). In the 
strong advection scheme, the strain term dominates near the corner of the 
cell, and the flame front cannot reach the corner. At d = 0.01, Figure [8^a) 
shows that incomplete combustion occurs near the corners of the cells, yet 
the flame front still manages to propagate. At d = 0.02, however, the flame 
front stops moving after t = 0.6. Note that if the level set stops moving, then 
the level set function forms a sharp layer. Here reinitialization is needed to 
alleviate the stiff level set function and keeps computation going. 

Figure [9^a) shows the graphs of s^ r (d) with A= 4, 6. In contrast to 
s^ s (d) > sl for any d > [15], sf r (d) decreases to zero when d is large 
enough. Figure Mb) shows the graphs of s s ^ r {A) with d = 0.01, 0.02. When 
A is small, (sl — dS) remains positive and s^ r is increasing. When A gets 
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Figure 9: (a) Plots of st = st(^) for strain G-equation with A = 4, 6. (b) 
Plots of s<r = st(A) for strain G-equation with d = 0.01, 0.02. 



larger, s?f r decreases and eventually drops down to zero. This agrees with 
the nonlinear phenomenon in turbulent combustion that high strain is the 
cause of flame quenching [U E] . 



5 Conclusion 

We have studied various G-equation models numerically, and evaluated the 
corresponding turbulent flame speeds in cellular flows. Based on the nu- 
merical results, we showed how the turbulent flame speeds are affected by 
viscosity, curvature or strain effect. Weak and strong bending effects of the 
speeds caused are observed in curvature and viscous G-equations. Quench- 
ing effect only appears in the strain G-equation. In future work, we plan 
to study turbulent flame speeds of G-equations in time dependent or three 
dimensional spatially periodic vortical flows. 



A Appendix: Surface Stretch Rate Formula 

In this appendix, we derive the surface stretch rate. A surface stretch rate 
formula in three dimensions is derived in [20j. Here we give an alternative 
formula in any dimensions and apply it in G-equation. 

Theorem A.l Suppose a smooth hypersurface in M d is moving in the ve- 
locity field V(x,t). Denote a the surface element area andn the unit normal 
vector of a point on the surface. Then the surface stretch rate is given by 

1 — = div(V)-n-DV-n. (A.l) 
a dt 
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Figure 10: An illustration of surface stretch in the proof of Appendix A. 



Proof: See figure 10 for the picture of the proof. Fix a time t and a point 
x on the surface, the surface can be locally approximated by its tangent 
plane. Let {m, . . . , rid—i} be an orthonormal basis of the tangent plane 
and ei, . . . , e^-i be infinitesimal scalars. Then the surface element can be 
presented by a rectangle whose sides are the vectors eim, . . . , ed-ind-i- The 
surface element area is 

a(t) = ei---e d _i. 

For 1 < k < d— 1, denote x k = x + e k n k the neighboring point of x of the 
rectangle. Then we say the rectangle is determined by the staring point x 
and neighboring points x%, . . . , Xd—%. 

After a time interval St, suppose the new locations of x, x k are x', x' k 
respectively. Then the surface element becomes a parallelogram determined 
by the staring point x' and neighboring points x[, . . . ,x' d l . Denote <5& = 
x' k — x', then the sides of the parallelogram are the vectors Si, ... , S^-i- The 
surface element area is 



a{t + 5 t ) = ^det(A T A), 

where A = [S\, ■ ■ ■ , Sd-i] is the matrix whose columns are the sides of the 
parallelogram. 

From now on we keep all calculations up to first order of St and omit 
higher order terms. The surface moves in velocity field V(x,t), then 

x = x + V(x, t)8 t , x' k = x k + V(xk, t)Sf 

=> S k = (xk -x) + DV ■ (x k - x)S t = (Id + S t DV) ■ (e k n k ). 
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Denote r] k = (I d + S t DV) ■ n k and N = [n\, ■ ■ ■ , n d _i], then 5 k = e^k and 



a(t + 5 t ) = e l ---e d ^dei(B T B), 
where B = [rj\, • • • , rja-i] = (Id + 5tDV)N. Then we have 

B T B = N T (I d + 5 t DV T ){l d + 5 t DV)N = I d ^ + 5 t N T (DV + DV T )N 
=> det(5 T J B) = 1 + 5 t tT(N T (DV + DV T )N) = 1 + 2<5 t tr(iV T WA0 

=> <j(i + <5 t ) = <r(i)(l + 5MN T DVN)). 
Hence the surface stretch rate is 

lim -^ a{t + 5t )- a ^ =tr(N T DVN)) 



S t ^o <r(t) 



St 



= n[DVni + ■■■ + n^DVnd-i- 
Note that {ni, . . . , n d _\, n} is an orthonormal basis of then 

njDVni + ■■■ + n T d __ x DVn d _ x + n T DVn = tr(DV) = div(V). 
We combine the last two equations and finish the proof. □ 



Remark A.l Result of 120^ in three dimensions reads: 

— — = (n • V)div(n) — curl(y x n) ■ n. 
a dt 

Indeed we can verify that (A.l) and (A. 2) are equivalent in IR 3 . 



(A.2) 



Corollary A.l Let V(x, t) be an incompressible flow and denote k = div(n) 
the curvature of the surface. If the surface moves in the velocity field V(x,t) 
and the normal direction with constant speed sl: 

dx 
~dt 

then the stretch rate is 

Ida _ 
a dt 



V(x,t) + s L n, 



-n ■ DV ■ n + slk. 



(A.3) 
(A.4) 



Proof: Substitute (A.3) into (A.l), then we have 
1 da 



a dt 



div(V) + Sidiv(n) — n ■ DV ■ n — slu ■ Dn ■ n. 



The first term is due to incompressibility of V. By some calculations, the 
last term is 0. □ 
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B Appendix: Iteration Scheme for Cell Problem 
of Viscous G-equation 

In this appendix, we prove the convergence of the iteration scheme for the 
cell (corrector) problem of viscous G-equation at large enough d: 

-ds L Au (k+1) + V(x)-Du (k+1) = - s L \P + Du {k) \ - V(x)-P, x G T n , 



#(*) = SL f \P + Dv,W\dx. 



(B.l) 



First we verify the solvability of (B.l). Denote L 2 er and Hp er the spaces 
of all mean zero and periodic functions in L 2 (T n ) and ff 1 (T n ) respectively. 
Since V(x) is assumed to be periodic, mean zero and divergence free, by 
Fredholm alternative theorem, the equation 

-Au + V{x) ■ Du = f , xGT™ 

has unique weak solution u G H^ er provided / G L 2 er . If G H, 



per 



per 



then the right hand side of (B.l) is in L 2 and there exists unique solution 



u 



£ Hp er for (B.l). Therefore given any G Hp er then we can 



construct a sequence {u^jfceN hi Hi, 



per " 



Theorem B.l The sequence {u^}ken in Hp er defined by the iteration scheme 
(B.l) converges provided d > ^/n/n. 



Proof: Replace the index k in (B.l) by k + 1 and take their difference: 
-ds L A(u^ - u^) + V(x) • £>( u (*+ 2 ) - «(*+*)) 
= (#(*+l) _ #(*)) _ SL Up + D U (H1) _ \ P + 
Multiply the equation by u^ fc+2 ^ — tt( fe+1 ) and take integration over T n : 



T™ L 



dx 

(fc+2) _ 



|p + _ |p + | ( u (*+2) _ u (*+l))dx. 



(B.2) 



Here we use the fact that V(x) is divergence free and u^ k+2 ^ — is mean 

zero. Recall the Poincare inequality: 



«||l2(T") < H^MIl^t™) > u ^ H^ 



7T 



per " 
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By Cauchy inequality, (B.2) implies that 

|2 



< || \P + W fe+1 ) I - \P + | || L2(Tn) \\u^ - n( fc+1 ) || i2(T „) 



||W fc + 2 ) - W fc+1 )||i 2(TIl) < ^||W fe+1 ) - W fc )|| L2(Tn) . 

If d > y/n/ir, then {Du^^jfcgN is contracting in L 2 (T"). By Poincare in- 
equality, {uW} fceN converges in Hp er . D 
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